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Abstract 

We compare the performance of several discretizations of simple 
pendulum equation in a series of numerical experiments. The stress 
is put on the long-time behaviour. We choose for the comparison 
numerical schemes which preserve the qualitative features of solutions 
(like periodicity). All these schemes are either symplectic maps or 
integrable (preserving the energy integral) maps, or both. We describe 
and explain systematic errors (produced by any method) in numerical 
computations of the period and the amplitude of oscillations. We 
propose a new numerical scheme which is a modification of the discrete 
gradient method. This discretization preserves (almost exactly) the 
period of small oscillations for any time step. 
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1 Introduction 



New but more and more important direction in the numerical analysis is 
geometric numerical integration [11], [12], [El [19] . Numerical methods within 
this approach are tailored for specific equations rather than for large general 
classes of equations. The aim is to preserve qualitative features, invariants 
and geometric properties of studied equations, e.g., integrals of motion, long- 
time behaviour and sometimes even trajectories (but it is difficult, sometimes 
even impossible, to preserve all properties by a single numerical scheme). 
"Although the apparent desirability of this practice might be obvious at first 
glance, it nonetheless calls for a justification" 

In this paper we perform a series of numerical experiments comparing the 
performance of several standard and geometric methods on the example of 
the simple pendulum equation. The equation itself is very well known but 
its discrete counterparts show many interesting and unexpected features, for 
instance the appearance of chaotic behaviour for large time steps [T] 
We focus our attention on the stabihty and time step dependence of the 
period and the amplitude for several discretizations of the simple pendulum 
(assuming that the time step is sufficiently small). We describe and explain 
small periodic oscillations of the period and of the amplitude around their 
average values. 

We confine our studies either to symplectic maps or to energy-preserving 
maps. It is well known that symplectic integrators are very stable as far as 
the conservation of the energy is concerned. Since the beginning of 1990s 
they are successfully used in the long time integration of the solar system 
[TTl ESI , see also [H H] . The reason is that using any symplectic scheme 
of nth order the error of the Hamiltonian for an exponentially long time is 
of the order 0{e"') where e is the constant step of the integration [S] [TT], [18]. 
Therefore, in studies of the long-time behaviour, symplectic algorithms have 
a great advantage at the very beginning. Fortunatelly, the class of symplectic 
integrators includes such well known and relatively simple numerical schemes 
as the standard leap-frog method and the implicit midpoint rule. In this 
paper we compare these classical methods with new geometric methods which 
preserve the energy integral. 

We also propose a new discretization (a modification of the discrete gra- 
dient method) which has some advantages: it is almost exact for small oscil- 
lations (even for large time steps) and keeps some outstanding properties of 
the discrete gradient method (e.g., its precision in describing motions in the 
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neighbourhood of the separatrix). 



2 Symplectic discretizations of Newton equa- 
tions 

We consider scalar autonomous Newton equations: 

^ = fi^) , (1) 
which can be written as the following first order system 

= P , P = f{^) ■ (2) 

The equations are integrable for any function / = f{ip) (in this case by 
integrability we mean the existence of the integral of motion, compare |26j). 
The energy conservation law reads 

+ V(,^) = E , /(^) = , (3) 

where E = const. The Hamiltonian is given by 

2 

H{p,q) = ^ + V{q) . (4) 

As an example to test quantitatively various numerical methods we will use 
the simple pendulum equation 

(p = —k sin if . (5) 

In this case the energy conservation law has the form 

-p^ — k cos ip = E . (6) 

The constant k is not important. It can be eliminated by a change of the 
variable t. In the sequel (in any numerical computations) we assume k = 1. 

By the discretization of we mean an e-family of difference equations (of 
the second order) which in the continuum limit e — > yields ([1]). The initial 
conditions should be discretized as well, i.e., we have to map ip{0) ipo, 
0(0) ^pq. 
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It is convenient to discretize ([2]) which automatically gives the discretiza- 
tion of p. Thus we have an e-dependent map p„) {ipn+i,Pn+i)- This 
map is called symplectic if for any n 

dipn+l A dpn+1 = difn A dpn ■ (7) 

The following lemmas give a convenient characterization of symplectic maps 
and we will apply them in the next sections. 

Lemma 1. The map {fn,Pn) ^ (v^n+i, Pn+i); implicitly defined by 

ipn+l - Vn = P{Pn, Pn+1, s) , Pn+1 - Pn = R{Vn, s) , (8) 

where P and R are differentiable functions, is symplectic if and only if 
dP OR _ dP dR 

dpn dipn dpn+l Oifn+l 

The proof is straightforward. Differentiating ^ we get 

difn+i - dipn = P,l dpn + P,2 dpn+1 , 
dpn+1 - dpn = R,l dipn + R,2 d^n+1 , 

(where the comma denotes partial differentiation). Then 
l + PaR^, , P,l+P,2 , 

= Z diPn + dpn , 

i — -r,2 -fx,2 i — -r,2 rL,2 

R,i +R,2 , 1 + P,i R,2 

dPn+l = p p d^n + -. p p dpn , 

provided that P,2 -R,2 7^ 1 (this condition means that the map defined by P, R 
is non-degenerate). Therefore 

difn+l A dpn+l = ^ p'^ p'^ d(fn A dpn ■ 

i — -r,2 ri,2 

Hence the map is symplectic if P,i R,i = P,2 P,2 7^ 1 which ends the proof. 
Lemma 2. The map {(pn,Pn) ^ (v?n+i, Pn+i), defined by 

(Pn+l - A{Lpn, e) + ipn^l = , Pn = fJ^oi^) '^n+l + B {ipn, s) , (10) 

is symplectic for any differentiable functions A, B. 
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In order to prove Lemma [2] we compute 

dpn+i = /iO dipn+2 + TB'dipn+l = 1^0 TA' d(pn+i - /io d(pn + TB' dipn+1 , 

where the prime denotes the differentiation and T denotes the shift. There- 
fore 

d(pn+l A dpn+1 = -yUo dlfin+l A dlfin ■ 

On the other hand d(pn A dpn = /io d{pn A d{pn+i, which ends the proof. 

3 Nonintegrable symplectic discretizations 

In this section we present some well known discretizations which preserve the 
symplectic structure of the Newton equations (compare [11], p. 189-190) but 
have no integrals of motion. 

3.1 Standard discretization 

The standard discretization of the simple pendulum equation 

= -ksm(pn (11) 

is non-integrable [26] . This discretization can be obtained by the application 
of either leap-frog (Stormer-Verlet) scheme or one of the symplectic splitting 
methods. It is interesting that we get the same discrete equation (ITT]) but a 
different dependence of Pn on ipn,^n+i (compare (ITBIl . (122]) ): 

Pn = h eke sm ipn , (12) 

where c = 0, |, 1. By virue of Lemma [2] standard discretizations are sym- 
plectic (for any c). 

3.2 Stormer-Verlet (leap-frog) scheme 

The numerical integration scheme 

Pn+l =Vn + \ef{iPn) , 

Vn+1 = ^n + ' (13) 

Pn+l = Pn+l + ^efi'fn+l) , 



is known as the Stormer-Verlet (or leap-frog) method (compare, e.g., 
Ehminating p^+i, we can easily formulate the Stormer-Verlet as a one-step 
method: 

(14) 

Pn+l =Pn + \e (/(V^n) + / (v^n + SPn + ^^^^(V^n))) • 

We can also formulate this method as 
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- = f{^n) , (15) 



Pn = ^ l^fKVn) ■ (16) 

In the simple pendulum case (/(</?) = —k sin y?) we recognize in the equations 
f fT5|) . (dnD the standard discretization ([n]), ([I2D with c = 1/2. 

3.3 Symplectic splitting methods 

The system (12]) belongs to the class of "partitioned systems" which have the 
form 

= giv^p) , p = Hv,p) > (17) 

where g, h are given functions of two variables. We can discretize such sys- 
tems in one of the following two ways: 

<fn+l = <fn + £g{y^n,Pn+l) , Pn+l = Pn + '^K'^n, Pn+l) , (18) 
<fn+l = <fn + egifn+lyPn) , Pn+l = Pn + eh{(pn+l,Pn) ■ (19) 

Both these discretizations are called either symplectic Euler methods [11] or 
symplectic splitting methods [2T]. In our case (see (E])) we have, respectively, 

V^n+l = ^n + ePn+1 , Pn+l = Pn + ^fi^n) , (20) 
(fn+l = (Pn + ePn , Pn+l = Pn + ^/(v^n+l) • (21) 

Finally, both ^ and ^ yield but instead of ^ we have 

<^n+l — V'n \ V^"+l ~ Vn /^^x 
Pn= £f{^n) or Pn = , (22) 



i.e., in the simple pendulum case we get f lT2]) with c = 1 and c = 0, respec- 
tively. 
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3.4 Implicit midpoint rule 



Any first order equation x = F{x) can be discretized using implicit midpoint 
rule (which coincides with the implicit 1-stage Gauss-Legendre-Runge-Kutta 
method, compare [H]). The first derivative is replaced by the difference 
quotient and the right hand side is evaluated at midpoint + Xn+i)- In 
the case of the simplest Hamiltonian systems, given by ([2]), we have: 

(fk+i = (pk + y{pk + Pk+i) , 

(23) 

In the special case of the simple pendulum we get 



^,+,-2^,+^,., _ (sin( '^''+f ^'' ) + sin C^'^+^^^-M ) , 

(24) 



£ 



2 ~ 2' 



Vk+l—fk _|_ Ic-Loi-r, Vk+l+'Pk 



Pk = \ + i^ek sin ^ 



The implicit midpoint rule has quite good properties: this is a symplectic, 
time-reversible method of order 2. The symplecticity follows directly from 
Lemma [H Indeed, fl23|l implies P,i = P,2 and -R,i = R,2- 



4 Projection methods 

Non-integrable discretizations can be modified so as to preserve the energy 
integral "by force", i.e., projecting the result of every step on the constant 
energy manifold. In principle, any one-step method can be converted into 
the corresponding projection method. In this paper we apply these proce- 
dures to the Stormer-Verlet (leap-frog) method. Therefore referring to the 
"standard projection" and "symmetric projection" we always mean standard 
(or symmetric) projection applied to the leap-frog scheme. 

4.1 Standard projection method 

There are given a first order equation x = F(a;), x G M^, any one-step numer- 
ical method = $e(a;„) (a discretization of the ODE), and a constraint 
g{x) = we would like to preserve. The standard projection consists in 
computing Xn+i '■= ^^(xn), and then orthogonally projecting Xn+i on the 
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manifold g{x) = 0, see [U]. This projection, denoted by Xn+i, yields the 
next step: x„ Xn+i- In other words, we define 

Xn+l = Xn+1 + \Vg{Xn+l) (25) 

where A is such that g{xn+i) = 0. 

Applying this approach to the simple pendulum ([5]) it is convenient to 
define 

X = (^V, = (<^, P) , (26) 

where uj = \/k. The above definition of p yields dimensionless components of 
X. li k = 1 (which is assumed throughout this paper), then this definition of 
p coincides with the previous one, see ([2]). The constraint g{x) = is given 
by dS]), i.e., 

g(x) = - cos V? - /i . (27) 
where h = E/uj"^. The equation fl25l) becomes 

iPn+l = + A sin ipn+1 , Pn+1 = (1 + ^jPn+l (28) 

and A is computed from 

^(1 + Xypl+i - cos{(fn + A sin (p^+i) = h . (29) 

In order to solve (1291) we use Newton's iteration Aj+i = Xj — AA^, where 

^ 1(1 + Xj)^pl^^ - cos(^n + \j sin^„+i) - h 
^ pI+i + sin^ 

and it is sufficient and convenient to choose Aq = 0. The approximated 
solution to fl29l) is given by A = limj^oo Aj. 



4.2 Symmetric projection method 

A one-step algorithm = $£(x„) is called symmetric (or time- reversible) 
if $_e = Equations of the classical mechanics are time- reversible, there- 
fore the preservation of this property is convenient and is expected to improve 
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numerical results. The symplectic splitting methods are not time-reversible 
while the Stormer-Verlet method and implicit midpoint rule are symmetric. 
The symmetry can be easily noticed in the form f|T3|) of the leap-frog method. 

The symmetric projection method preserves the time-reversibility. The 
method is applied under similar assumptions as standard projection (addi- 
tionally we demand the time-reversibility of and consists of the following 
steps |21[in]: 

Xn = Xn + X'^giXn) , 



Xn+1 — Xn+1 + AV5'(x„_|_i) , 

where we assume g{xn) = and compute the parameter A from the condition 
gixn+i) = 0. 

5 Integrable discretizations 

Throughout this paper by integrability we mean the existence of an integral 
of motion. The Newton equation ([T]) has the energy integral Its dis- 
cretization is called integrable when it has an integral of motion as well. In 
the continuum limit this integral becomes the energy integral, so it may be 
treated discrete analogue of the energy. 

5.1 Standard-like discretizations 

Standard-like discretizations are defined by [26] 



Pn+1 =Pn+ £F{ipn, e) 

where F has to satisfy F{ipn, 0) = f{fn)- For a given / there exist inifinitely 
many functions F satisfying this conditions. All of them are symplectic, 
which can be easily seen applying Lemma [1] with P,i = R,2 = 0. Similarly as 
in Section [3] we obtain from (132|) : 



Xn+1 = ^e{Xn) 



(31) 



(32) 



(33) 



Pn = 



e 



eF{^n,e) 



e 
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We are interested in integrable cases, i.e., in discretizations preserving 
the energy integral. Suris found that two standard-hke discretization of the 
simple pendulum are integrable [25], [26] : 



^ ( /ce^sinoj^ \ 
LPn+i - 2v9„ + v^n-i = "2 arctau , (34) 

\ Z ~r COS / 

ifn+i - 2v5„ + = -4 arctan . (35) 

\4 + te^ COSipnJ 

The equation ( 134|) . referred to as the Surisl scheme, has the integral of motion 
given by 

1 / 2 sin '^"+^~'^" V 1 

Ei = -\ — ^ J - -k (cos (pn + COS (fn+l) , (36) 

or, in terms of (fn and Pn, 



1-cosepn I, , , . 

El = -k (cos V^n + COS(v3n " £Pn>) ■ (37) 



The equation ( ]35l) . referred to as the Suris2 scheme, has the following integral 
of motion 

^2 = - I I - A; COS ^ , (38) 

which can be expressed in terms of and p„ as follows 
— (^1 - cos — j - fc cos(v3„ ^ 



E2 = — [I — cos ) — ^ cos(v9„ ^) . (39) 



One can verify the preservation of these integrals by direct computation. 



5.2 Discrete gradient method 

The discrete gradient method [211 1221 [23] is a general and very powerful 
method to generate numerical schemes preserving any number of integrals of 
motion and some other properties [20]|. However, this method in general is 
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not symplectic. In this paper we need to preserve one integral (the energy) 
and the system is hamiltonian, compare (jl]), 

dH dH 

In such case the discrete gradient method reduces to the following simple 
scheme. Left hand sides of the formulas (j^Ol) are discretized in the simplest 
way (difference quotients) while the right hand sides are replaced by the so 
called discrete (or average) gradients: 

e ~ Aif ' e ~ Ap - ^ ^ 

The discrete gradient VH = (^^, of a differentiable function H{ip,p) 

by definition (see [21]) satisfies the condition 

AH AH 

The explicit form of VH is, in general, not unique. One of the possibilities 
is the coordinate increment discrete gradient [T5] 



AH H{ipn+l,Pn) - H{ipri,Pn) Aif if(v?„+l, p„+l) - if (</?„+!, Pr, 



A(p ifn+l - Ap Pn+l - Pn 

(43) 



Other possibilities are, for instance, mean value discrete gradient |2T] and 
midpoint discrete gradient [8]. All these definitions coincide in the case 
H{ip,p) = T{p) + V{ip). In such case VH = VT + VV^, where 

Pn+l - Pn ' ^n+1 " V^n 

Thus we have got the discrete gradient scheme: 

Pn+l + Pn " ^n 



2 e 

Pn+l - Pn ^ Vj^n+l) - Vj^n) 
£ y^n+1 - <fn 



(45) 
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This numerical scheme can be also obtained as a special case of the modified 
midpoint rule [TB] . The system (HSj) can be rewritten as the following second 
order equation for (^9^ plus the defining equation for 

V?„+l - 2ifn + V^n-l _ 1 [ V{ifn+l) " V {ifn) _^ V{ipn) - V{(pn-l) 



(46) 

Pn = h -e 



Substituting V{ip) = —k cos ip we get the simple pendulum case. Multiplying 
both equations ( l45l) side by side, we easily prove that the system ( l46l) has 
the first integral 

E=IpI + V{if^) (47) 

which exactly coincides with the hamiltonian (jlj) evaluated at ifn, Pn- Note 
that the integrals of motion fl37|) . fl39l) coincide with (jl]) (where V{ip) = 
—k cos if) only approximately, in the limit e —>■ 0. 



6 A correction which preserves the period of 
small oscillations 

The classical harmonic oscillator equation (p + u'^if = admits the exact 
discretization ([H], compare also [HEl]), i-e., a discretization such that the 
solution ip{t) evaluated at ne equals ipn (for any e, and any n): 

(fn+l - 2ipn COS eUJ + (fin-l = 0, 

. . . (48) 

Pn = - COS Ue) . 

smue 

The energy is also exactly preserved, i.e., 

E = IpI + l^Vn (49) 

does not depend on n (which can be easily checked by direct calculation). The 
existence of the exact discretization of the harmonic oscillator equation has 
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been recently used to discretize the Kepler problem (preserving all integrals 
of motion and trajectories) [5]. 

We consider the class of Newton equations ([1]). Let us confine ourselves 
to equations which have a stable equilibrium at ip = 0, i.e., /'(O) < 0. Then 
V = V{ip) has a local minimum at v? = 0, i.e., ^'(0) = /(O) = 0. We denote 



^0 = VV^"(0) . (50) 

Thus 

V{y,) = Vo + ^^oV + • • • , (51) 

and small oscillations around the equilibrium can be approximated by the 
classical harmonic oscillator equation with u = ujq. 

Do exist discretizations which in the limit y^n ~ (e is fixed) become 
exact? Known discretizations, including those presented in this paper, do not 
have this property. Fortunatelly, we found such discretization by modifying 
the discrete gradient method. It is sufficient to replace e by some function 
6 = 6{e) in the formulae psl) . The form of this function will be obtained by 
the comparison with the harmonic oscillator equation (in the limit ^ 0). 

We linearize the equations fH6|) (with e replaced by 6) around = 
(i.e., we take into account fIST]) ). Thus we get 



(52) 

Vn = ^ h -UJQd (V9„+1 + ipn) , 

which is equivalent to 



(53) 



We compare ( HHl) with ( l53i) . Both systems coincide if and only if 

^ = cos ecu , — ^ = . (54) 
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Solving the system we get 




) 



(55) 



Therefore, we propose the following new discretization of the Newton equa- 
tion ([1]), ([3]) (modified discrete gradient scheme): 



where 5 is defined by ( 1551) (and ojq is given by (!50|) ). This discretization 
becomes exact for small oscillations for any fixed e. It means that for ipn^ ^ 
the period and the amplitude of the approximated solution should be very 
close to the exact values (even for large e\). In the next sections we will verify 
this point experimentally. 

7 Numerical experiments 

We performed a number of numerical experiments applying the numerical 
schemes presented above. The initial data were parameterized by the velocity 
Po while the initial position was always the same: y^o = 0. In the continuous 
case ([5]) we have 3 possibilities: oscillating motion (|po| < 2), rotating motion 
(IpoI > 2) and the motion along the separatrix (po = ±2), from = to 
(asymptotically) (p = ±7r. The (theoretical) amplitude Ath for the oscillating 
motions can be easily computed from the energy conservation law (jH]) (where 
k = 1, i.e., = - cos Ath): 



In particular, we performed many numerical computations for the following 
initial data: 



• Po = 0.1, then Ath ~ 0.03184437r ^ 0.1000417 (small amplitude) 

• Po = 1.8, then Ath ~ 0.7128 677r ^ 2.239539 (very large amplitude). 



(fn+l - 2(pn + V?n-1 

6^ 




(56) 





(57) 
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To estimate the actual amplitude of a given discrete simulation we apply 
the following procedure: if ip^ is a local maximum of the discrete trajectory 
(i.e., ifm > ^Pm-i and > Pm+i), then we estimate the maximum of the 
approximated function by the maximum of the parabola best fitted to the 
following five points: (pm-2, ^m-i, ^m, ^m+i, The analogical procedure 

is done also at local minima (we take the absolute value of the obtained 
minimum). Thus we obtain a sequence of the amplitudes, Ajy. The index 
N is common for all extrema (maxima and minima), and on some figures 
we denote it by N1/2 (the number of half periods) to discern it from N (the 
number of periods). 

Every numerical scheme used in the present paper yields a discrete tra- 
jectory with rather stable amplitude. It is not constant but oscillates in a 
regular way around an average value: 



where both the average amplitude A and relative (dimensionless) oscillations 
Oat can depend on the time step e and on the initial velocity po, i.e., A — 
A{po,e) and ajv = Q;jv(po,£)- Of course, both A and ajv differ for different 
numerical schemes. 

In a similar way we estimated the period of discrete motions. The exact 
periodicity {(pk+n — for some n) is a rare phenomenon and, of course, 
we did not observe it. To define the approximate period we fit a continuous 
curve to the discrete graph, estimate zeros of this function, and compute the 
distance between the neighbouring zeros. 

Suppose that iprn^m+i < for some m. It means that one of the ze- 
ros, say zjsi, lays between c/?^ and (pra+i- We estimate it by zero of the 
interpolating cubic polynomial based on the points ^Pm-ii fm, fm+i, 'fim+2 
(another natural, but less accurate, possibility could be a line joining and 
Pm+i)- Then, denoting subsequent estimated zeros by zn {N = 1,2,3,...) 
and zq = Po = 0, we define 

Tn — Z2N — Z2N-2 , (59) 

which we take as an estimate of the period. 

Our numerical experiments have shown that Tn is not exactly constant 
but oscillates with a relatively small amplitude. The average value of is 
constant with high accuracy (see the next section). Therefore we have 



An = A{1 + aN) , 



(58) 



Tn - T{1 + Tn) , 



(60) 
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where both the average period T and relative (dimensionless) oscillations 
Tjv can depend on the time step e and on the initial velocity po, i.e., T = 
T{pQ,e) and = tn{po,£)- Moreover, T and tn essentially depend on the 
discretization (numerical scheme). 

The amplitude of small oscillations is defined in a natural way 

T{e,po) := max |rAr(£:,po)| , Oi{e,po) := max\aN{e,po) . (61) 

Fortunatelly, |r7v| and \ai\f \ oscillate (as functions of N), with small ampli- 
tudes, in a very regular way. Thus we can estimate r(e,po) and a{e,po) 
considering a series of, say, 40 local extrema of tn and on, and taking an 
average value. 



8 Periodicity and stability 

Discrete trajectories generated by symplectic or integrable schemes consid- 
ered in our paper are stable for e which are not too large (for very large e one 
can observe chaotic behaviour, [3 [29]). We confine ourselves to sufficiently 
small e, i.e. e ^ 0.5, but sometimes (for pq < 1.5) we can take even e ^ 1. In 
this region the motion is very stable and both the average period T and the 
average amplitude A are well defined. The average amplitude is computed 
simply as 

^ M-l 

A,„,(iV,M) = -5^|A^+,| , (62) 

j=0 

where we usually assume M = 50. The definition of the average period is 
similar. In many cases we use the formula 

Tavg{N, M) = ^ {zn+2M - Zn) , (63) 

where the dependence (very essential!) on e and po is omitted for the sake 
of brevity. Note that = Tavgi^N — 2, 1). Computing Tavg it is necessary 
to choose M arbitrarily, we usually take M = 20. Sometimes we denote 
N = Nq to point out that the average is taken over indices greater than Nq. 

Considering very long discrete evolutions (many thousands of periods) we 
use another definition of the average period. Namely, we average Tavg{N, M) 
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over some range of the parameter M {K < M ^ L): 

1 ^ 

Ta.g{N,K,L) = j—j^ J2 T,.g{N,M). (64) 

M=K+1 

Usually we assume K = 100, L = 200. 

All discretizations considered in the present paper are characterized by 
very high stability of the period and the amplitude. One can hardly notice 
any dependence of Tavg and Aavg on A^, even when testing very large (like 
10^, 10^ or 10^), and Tavg is even more stable. 

As a typical example we present long-time behaviour of the Surisl scheme, 
see Fig. [T] and Fig. [21 where we used the definition 0631) with M = 20. An 
interesting phenomenon is associated with changing M. The pictures for 
different M usually are very similar but the amplitude of oscillations becomes 
smaller and smaller for larger M (compare Fig. [3|, where M = 1, with Fig. [2l 
where M = 20). 

Table [T] shows how stable are periods of the oscillations. Maximal T^v 
is defined as maxj+ioo<Ar^j+2oo for either J = or J = 1.8 ■ 10^. Min- 
imal values and the average are taken over the same range of values. The 
standard error of the average is about 5.7 • 10~^ (the maximal error is about 
10^^). Therefore, the average period is practically constant for all studied 
discretizations. The Surisl scheme is exceptionally stable. In this case any 
variations of the period are well within the error limits and we did not ob- 
serve any dependence of Tavg{N, M) on A^. Taking into account the observed 
stability of the period, throughout this paper we identify the average period 
with T = T,,g (0,20). 

The observed stability of the period (for symplectic and integrable dis- 
cretizations) is in sharp contrast with the results given by standard (non- 
symplectic and non-integrable) numerical methods. For instance, the most 
popular (explicit) 4th order Runge-Kutta scheme yields the period noticeably 
decreasing in time (see Fig. H]). For small A"o we get reasonably good estima- 
tion of the period (interpolating the discrete curve we get T = 11, 64602 for 
A'o = 0, which is quite close to the theoretical value Tth = 11, 65758528. From 
among our discretizations only both gradient schemes produce comparable 
(even a little bit better) results, namely the discrete gradient scheme yields 
T = 11.64698. However, for larger A^o the Runge-Kutta method yields worse 
and worse estimation of the period (in fact this is an exponential decrease, 
although very slow) while both gradient methods remain stable for very long 
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time, compare Tabled! In this particular case (po = 1-95, e = 0.2) the error 
produced by the Runge-Kutta method becomes greater than the errors of all 
methods considered in this paper beginning from A^o ~ 2000. 

Numerical experiments show that the oscillations of the period and the 
amplitude are very small. For e — * we have r(£,po) 0, up to the round- 
off error. The largest values of T{e,po), obtained for both projection methods 
(for large e and small po), are of order 0.2. All other discretizations yield 
oscillations smaller by one or two orders of magnitude (even for large e). A 
typical picture is given at Fig. [5] representing r(e,po) for po = 1.8. 

9 Why the period and the amplitude oscillate 
in a very regular way? 

In a large range of parameters the oscillations tn are very regular and their 
amplitude is greater than numerical errors by several orders of magnitude. 
This phenomenon turns out to be caused mainly by systematic numerical 
by-effects. 

Our explanation is associated with the above procedure of estimating ze- 
ros. In general, the period T = Tavg and e are incommensurable. Therefore 
the relative position of zn between and ipm+i depends on N. We con- 
jecture that the periodic phenomena one observes at Fig. El Fig. [TJ Fig. [HI 
Fig. [9], Fig. [To] and Fig. [11] are associated with the properties of the real 
number T/e, namely, with the approximation of T/e and T/{2e) by rational 
numbers. 

We begin with a simple definitions. Given T, e G M (T > £ > 0) and 
K E N we define: 

KT KT 
■= Mk , i^K ■= — Lk , (65) 

such that -0.5 < fix ^ 0.5, -0.5 < uk ^ 0.5 and Mk,Lk G N. In other 
words, for a given K we take Mk such that Mk/K is the best rational 
approximation (with a given denominator K) of the real number T/e, and 
Lk/ K is the best rational approximation (with the denominator K) of T/e. 
For given T,e,K the formulas ([65!) define uniquely fix, ^k, Mk, Lk- The 
following lemma can be derived directly from the above definitions. 

Lemma 3. Suppose that T > e > are given. 
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1. If \jjK + I^jI < 0.5, then Mk+j = + Mj and Hk+j = ^-k + /^j- 

2. If \i^K + < 0.5, then Lk+j = Lk + Lj and uk+j = i^k + ^j- 

3. If \i^k\ < 0.25, then Mk = 2Lk and hk = 2uk- 
4- If K is even, then Mk/2 = Lk and fiK/2 = ^k- 

Corollary 1. If uk ~ 0, then fix ~ and, for K even, also fiK/2 ~ 0. 

If fiK ~ 0, then the configuration of zjy, ip^, ipm+i practically repeats after 
every K periods. Therefore it is natural to expect some periodic recurrences 
with the period KT. In particular, t^+k ~ tn for any A^. 

To obtain a "good" approximation we usually demand at least fiK < 
0.01. Sometimes, especially for small K (e.g., K ^ 5), interesting effects can 
be observed also for larger (but, anyway, fix < 0.1): the graph of the 
function N apparently splits into K "discrete curves" (Tn and Tm 

belong to the same curve if = M (modi^')). 

Similar considerations can be made for the oscillations aN of the ampli- 
tude. In this case the period is T/2 and "good" approximations correspond 
to uk ~ 0. 

Example 1 (leap-frog scheme, e = 0.05, po = 1.8, T ^ 9.1254146). We 
compute T/e ^ 182.508291 and easily check that fi2 ~ 0.017, fi^g ^ —0.011, 
/i6i ^ 0.0058, /ii2o ~ -0.0051, /ii8i ^ 0.00067. Fig. \E confirms that the 
characteristic "time scales" responsible for the pattern of the oscillations are 
2, 120, and 181, indeed. 

The period 2 corresponds to oscillations between two sinusoid-like curves. 
Namely, Tjv belong to the first "sinusoid" for N odd, and to second "sinusoid" 
for N even. Both discrete curves are periodic with the period 120. Actually, 
the whole picture seems to have the translational symmetry with the period 
60. The difference between T^+m and T^ is quite large (in this sense 60 is 
not a period, indeed), however T^ lays between T/v+59 and T/v+gi- 

The next period, 181, is more dificult to be noticed and corresponds to 
more subtle effects, like the configuration of points near intersections of both 
"sinusoids" which approximately repeats every three " sinusoid" -half-periods. 

Similarly, we compute z/4 ^ 0.017, 1/59 ~ —0.0054, z^igi ^ 0.00034 and 
^240 ~ —0.0051. On Fig. \2\ we recognize four discrete curves, periodic with 
the period 24O. The whole picture has the period 60 but looking closely on 
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some details (e.g., at peaks or at intersections) we can also notice another 
periodicity with the period 181. 

Finally, we point out that all equalities suggested by Lemma\^hold (e.g., 

/i61 = yU2 + /i59, ^^240 = ^^59 + ^^181, A^SQ = 2Z/59, /i4 = V2 etc). 

Example 2 (leap-frog scheme, e = 0.1, po = 0.05, T ^ 6.28155042). T/e ^ 
62.815504 and we check that /ig = 0.078, /in = -0.029, IJ27 = 0.019, /igs = 
—0.011, fiQ5 = 0.0078, yUios = —0.0031. Fig. [3 does not look so regularly 
as Fig. O Note that fix o,re now relatively large, the first fix smaller that 
0.01 has the index K = 65 and the next one is K = 103. However, a closer 
inspection reveals similar features in both figures. We have five sinusoid-like 
curves (periodic with the period 65). The distance between them is 13 but the 
difference between T/v+13 and T/v is large. Note that the period 103 f« 8 x 13, 
so points of only every eighth '^sinusoid" practically coincide. 

The other periods (K = 11,27, 38j can be derived from 103 and 65, 
namely: 38 = 103 — 65, 27 = 65 — 38, 11 = 38 — 27. They can be noticed 
on Fig. [2| as well. For instance, the lowest points (Tjsi between 6.28155037 
and 6.28155038) have N = 6, 17, 22, 33, 44, 49, 60, 71, 82, 87, 98, the distances 
between them are given by AN = 11,5,11,11,5,11,11,11,5,11 (note that 
11 + 11 + 5 = 27;. 

To explain regularities on Fig. we compute z/5 = 0.039, z/22 = —0.029, 
1/27 = 0.00 93, 1/49 = -0.020, z/76 = -0.011, u^s = -0.0015, z/130 = 0.0078 
and also i/645 = 0.000 1 0. In this case the structure is also quite compli- 
cated because we have several candidates for periods. Some of them admit a 
clear interpretation. Joining every fifth point we get five sinusoidal curves 
with the period 130. Thus the distance between neighbouring "sinusoids" is 
26 which is very close to the period 27. The subsequent minima are at N = 
3, 25, 52, 79, 106, 128, 155, 182, 209, therefore AN = 22, 27, 27, 27, 22, 27, 27, 27 
(note that |z/22| is also relatively small). Looking at configurations of points 
near every minimum we can notice a distinct periodicity with the period 103. 

Example 3 (Surisl scheme, e = 0.1, po = 0.05, T ^ 6.29723795). In this 
case the structure of Fig. [JO is extremaly simple (a single discrete curve). It 
can be explained by the non-existence of any "small" periods. The smallest 
one, distinctly seen at Fig. [721 is 36. Namely, /ise = 0.0057, fi^^ = 0.0050, 
/^i8i = 0.00069. The period 181 is even more exact than the period 36 (fiisi 
is much smaller than fise). Therefore after every five basic periods (181 ~ 
5 X 36; the periodicity improves. 
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Fig. [71] consists of two intersecting discrete curves (periodic with the pe- 
riod 72), because z/2 = —0.028 is relatively small and z/72 = 0.0057. Actually 
the important point is that z/3 = 0.46 is much greater than |z/2|. Note that 
fi2 = —0.055 is also not very large but /i3 = —0.083 is of the same order. 
The whole structure has the period 36 but ( similarly as in Example [2]j the 
difference between Ajv+sg and is quite large, Ajq is close to A^^^.^ O'^'^d 
^^^35 = 0.017, z/37 = —0.011/ Moreover, we have the period 181, quite 
accurate (vi^i = 0.00034j. This periodicity can be noticed by looking at the 
minima or at points where the discrete curves "intersect" . 

Similar remarks concern the case presented at Fig. [H Fig. [21 Fig. [31 where 
T ^ 11,88884005 and /xg = -0.0044, /i448 = 0.0035, ^457 = -0.00093. The 
patttern on any of these figures consists of nine discrete curves and is periodic 
with the period close to 457. 

The behaviour described on the above examples is typical and similar 
periodic phenomena can be observed for other discretizations and for other 
choices of parameters except very small values of e (e.g., e ^ 0.01) when 
periodic oscillations are comparable or smaller than the round-off error (then 
the oscillations become chaotic with a very small amplitude). 

10 Numerical estimates of the amplitude and 
the period 

All discretizations considered in this paper are characterized by very good 
stabihty of their trajectories. Therefore, such quantities as average period 
and (in the case of oscillating motions) average amplitude are well defined 
for every discretization (provided that e is not too large, it is sufficient to 
assume e ^ 0.5). 

10.1 Average amplitude 

Relative errors for the average amplitude are presented in Table [2] (for e = 
0.02 and e = 0.5). They were computed as differences between the numerical 
results and exact amplitudes given in terms of elliptic functions. One can 
immediately see that in any case the best results are given by both gradient 
schemes (and the worst ones are given by Surisl and Suris2 schemes). The 
relative error of the leap-frog and Suris' methods practically does not depend 
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on pq. The accuracy of gradient methods increases for larger pq, both for e = 
0.02 and e = 0.5. For small e (e.g., e = 0.02) also projection schemes yield 
very small errors, like 10~^ or 10~^ (similar as gradient methods). However, 
for some po their accuracy is very high (e.g, for po = 1.6) while for some 
other Po ^ relatively worse (e.g., for po = 0.8). 

The implicit midpoint rule is comparable to gradient methods but only 
for small po (e.g., po < 0.1). The leap-frog method, both Suris' discretization 
and (for po > 1.6) the implicit midpoint rule yield much larger errors (by 4 
orders of magnitude). 

For greater e (e.g., e = 0.5) the differences between the studied methods 
are much smaller (they differ at most by 2 orders of magnitude). Gradient 
methods are most accurate. The implicit midpoint rule has similar accuracy 
for Po < 1.2 while projection methods are not much worse for po > 1.8. 
Leap-frog method and both Suris' methods have larger relative errors for 
any pq. We point out, however, that even those "large" errors are not so bad 
(only several percent) with the exception of po approaching 2 (when these 
discretizations fail to reproduce properly even the qualitative behaviour). 

Fig. [IS] illustrates the dependence of the average amplitude on e for po = 
1.8. Gradient methods and (especially for e < 0.3) projection methods are 
most accurate. 

10.2 Average period 

Relative errors for the average period are presented in Table [3] and also in 
Table H] (in both cases for e = 0.02 and e = 0.5). For po < 0.5 all dis- 
cretizations except the modified discrete gradient method have similar rela- 
tive errors (Surisl scheme is the worst among them). The modified discrete 
gradient methods is much better (for po ~ its error is smaller by 4 orders 
of magnitude, at least), compare Fig. [12] (po = 0.1) and Fig. [H] (po = 0.02). 

Then, with increasing po, all discretizations become to have similar ac- 
curacy with two very interesting exceptions: leap-frog and implicit midpoint 
schemes have a kind of "resonance values" for which their accuracy is much 
better than the accuracy of all other method. Fig. [TJ] shows how accurate is 
the leap-frog scheme for po = 1.21 and for practically any e. There are shown 
also next two discretizations: implicit midpoint and modified discrete gradi- 
ent, much worse (for this value of po) than leap-frog (other discretizations are 
even less accurate). Implicit midpoint scheme has an analogical "resonance 
value", namely po ~ 1.6. It is worthwhile to point out that, surprisingly. 
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projections applied to the leap frog method have strong negative effect on 
the accuracy of the average period for 0.8 < po < 1.8, especially for larger e 
(e.g., £ = 0.5). 

If Po approaches 2, then both gradient methods become more accurate 
than other methods (only for small e the projection methods are better). 
For pq very close to this limiting value the accuracy of all methods decreases 
rapidly, and the leap-frog method and both Suris' methods produce rotating 
motions instead of oscillations, see Table HI The closest neighbourhood of 
the separatrix {po = 2) is discussed in more detail below. Here we remark 
only that, for po slightly greater than 2, the implicit midpoint method fails to 
reproduce rotations and has wrong qualitative behaviour (i.e., oscillations). 

In the case of rotating motions the relative error of the average period is 
very similar for all considered methods except the discrete gradient scheme 
which is better by one or two orders of magnitude. 

11 Interesting special cases 

In this section we briefly present several points which seem to be encouraging 
to further studies. 

11.1 Extrapolation e ^ 

For all studied discretizations we expect 

\imT{£,po) = TthiPo) , \imA{e,po) = A/*(Po) (66) 

where Tthipo), Ath{po) do not depend on the discretization and are equal to 
theoretical values computed from the analytic formula (in terms of elliptic 
functions), compare Fig. [T2l and Fig. [T3l 

Let us analyse quantitatively the case presented at Fig. [T2] (the exact 
period is Tth ~ 6.28711783). Fitting 3rd-order polynomials (very close to 
parabolas, in fact) to twelve points {e = 0.01, 0.02, . . . , 0.11, 0.12) we get 

T = -0, 03867e3 + 1, 310512e2 _ 0, OOOlOSOe + 6, 28711875 (Surisl) , 

T = -0, 00909e3 + 0, 5240536^ - 0, 0000247e + 6, 28711805 (Suris2) , (67) 

T = -0, 00475e3 - 0, 2602426^ _ 0, 0000130e + 6, 28711794 (leap-frog) . 

The last terms estimate the exact period quite well. Taking 10~'' as a unit 
we compute their absolute errors as: 9.2, 2.2 and 0.9, respectively. They 
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are comparable with the errors at e = 0.001 (given by 13.1, 5.2 and —2.6, 
respectively). The errors at e = 0.01 (namely, 1307.2, 523.3 and 260.6) are 
higher by two orders of magnitude. The modified discrete gradient scheme 
(with the (5-correction) beats all other discretizations: its error at e = 0.01 is 
only 1.3 (in the same units). 

11.2 The neighbourhood of the separatrix 

The separatrix is a border between oscillating and rotational motions. Ta- 
ble H] presents the values of the period for motions near the separatrix, i.e., 
Po ~ 2. This is certainly the range of parameters most difficult for accurate 
numerical simulations. The gradient schemes and projection methods yield 
satisfying results, especially for small e, and are much better than all other 
methods. For rotating motions very close to the separatrix even projection 
methods (especially the symmetric projection) become less accurate and only 
gradient methods yield relatively good quantitative results, see Table HI 

The other discretizations produce wrong results (in the neighbourhood 
of the separatrix) even qualitatively. Namely, the leap-frog and both Suris' 
schemes begin to simulate rotating motions for po < 2 (e.g., for po = 1.99 if 
e = 0.5, and for po = 1.99999 if e = 0.02), while the implicit midpoint rule 
produces oscillating motions for Po > 2 (e.g., for pq = 2.000001 if e = 0.02, 
and for po = 2.001 if e = 0.5). Even in the case of good qualitative behaviour 
these methods yield very large relative errors, especially for larger e (for 
e = 0.5 and |po — 2| ^ 0.001 leap-frog, imphcit midpoint and both Suris' 
schemes yield relative errors like 30% — 70% and more. 

If Po = 2, then (in the continuous case) we have the motion along the 
separatrix, i.e., ^ vr for t —>■ oo. For larger e (e.g., e = 0.2) this behaviour 
is not reproduced by any discretization. Interesting results are given by both 
gradient schemes, see Fig. [T7| {e = 0.2). The standard gradient scheme pro- 
duces oscillations, but after three periods one rotation is performed. The 
modified discrete gradient scheme gives a strange motion: first oscillations 
(two periods), then backward rotation (3 periods), forward rotation and the 
return to oscillations. This picture depends on e and the round-off error 
chosen. In any case, for both gradient schemes, we have a number of chaotic- 
looking switches between oscillations and rotations in both directions. Qual- 
itatively this behaviour may be considered as satisfying. It reflects the fact 
that the equilibrium at = vr is unstable. In the same time, the projective 
discretizations (quite good at qualitative description of motions near the sep- 
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aratrix) produce relatively slow rotational motion (similarly as the standard 
leap-frog method and both Suris schemes). However, for very small e (e.g., 
e ^ 0.00025) the symmetric projection method seems to have the proper 
qualitative behaviour and is much better than other considered numerical 
schemes, see Fig. [181 

11.3 Advantages of the new method 

The discrete gradient method with (5-correction turned out to be very effi- 
cient as far as the numerical estimation of the period (for relatively small 
amplitudes) is concerned. The range of these "small" amplitudes is quite 
large, up to ~ 7r/4, which corresponds to po < 0.8. Thus it contains also 
the cases which cannot be approximated by the linear oscillator. Even for 
Po ~ 0.8 the new method is several times better than the best of other con- 
sidered schemes, and for smaller po it becomes better even by 4 orders of 
magnitude (e.g., for = 0.02 the errors of other discretizations are greater 
by the factor at least 0.5 ■ lO'', see Table [3]). 

Fig. [12] (po = 0.1) shows how precise is the period given by our new 
method in comparison to the period given by other numerical schemes. Sim- 
ilarly, Fig [TH presents the relative error for po = 0.02 and a large range of e. 
We see that even for e = 1 the relative error is only 10^^! For small e the 
error is 10~^ and less. 

Our method works very well also for larger amplitudes, but for po larger 
than 1,4 the discrete gradient method is better, and the leap-frog scheme 
and implicit midpoint are unbeatable around their "resonance" amplitudes 
(po ~ 1, 2 and po ~ 1, 6, respectively). In the case po > 2 the delta correction 
have negative influence on the accuracy of the gradient discretization (which 
is the best for rotating motions). However, the accuracy of the modified 
discrete gradient method is on the same level as the accuracy of all other 
considered methods. 

In the close neighbourhood of the separatrix the modified discrete gradi- 
ent scheme behaves similarly to the discrete gradient method and its qualita- 
tive behaviour is perfect. What is more, also the quantitative results are very 
good (compare Table [1]). Fig. [16] compares the behaviour of our method with 
the leap-frog and implicit midpoint schemes for po = 2.000001. The points 
generated by the modified discrete gradient method practically coincide with 
the exact solution (the relative error of the period is 0.59%), almost as good 
result as that given by the discrete gradient scheme (the error is 0.25%). 
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The leap-frog scheme produces good quahtative behaviour but with the pe- 
riod two times smaller than the exact one. The implicit midpoint scheme 
gives wrong qualitative result: oscillations instead of rotation. 

12 Conclusions 

All methods considered in this paper are characterized by very high stability 
of periodic motions they generate (provided that e is not too large). The 
average period is practically constant (with the accuracy close to 10"'' or 
better) for a very long time (we checked even several millions of periods). 
The period and the amplitude perform regular small oscillations (they are 
relatively larger for both projection methods). The periodic character of 
these oscillations turns out to be of a systematic origin and we explained it 
considering rational approximations (with possibly small denominators) of 
the real number T /e. 

The main aim of this paper was the comparison of several numerical 
schemes. The standard leap-frog method, although non-integrable, is quite 
good when compared with typical integrable discretizations. Its performance 
should be enhanced by use of projection methods which impose the conserva- 
tion of the energy integral. The projections work very well for small values of 
the time step (e.g., the symmetric projection gives excellent results simulat- 
ing the motion along the separatrix), while for larger time steps they produce 
relatively large fluctuations of the period and the amplitude. In any case the 
projections produce much more accurate values of the average amplitude. 
The average period is of the same order, or even worse (in comparison to the 
standard leap-frog method). 

Surprising resonances occur for po ~ 1-21 (for the leap-frog method) 
and po = 1-6 (for the implicit midpoint rule). In the neighbourhood of 
these "resonance" values these methods have exclusively high accuracy of the 
estimated period (practically for any e), much better than all other methods. 
It would be interesting to explain this phenomenon. 

Discretizations found by Suris |i2^ are very stable but, in the same time, 
they have relatively large errors as compared to other numerical schemes. 
This is surprising because these methods are both integrable and symplectic. 
In this case the error (i.e., deviation from the exact solution) seems to be of 
a systematic origin. We plan to construct appropriate modifications of Suris' 
discretizations in order to enhance their precision without destroying their 
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stability. 

The discrete gradient method is (for any e and any po) among the most 
accurate methods. For rotating motions this is certainly the best method. 
We proposed a modification of the discrete gradient method which proved to 
be quite successful, especially when applied to simulate oscillating motions. 
Our new method is extremaly efficient for small oscillations. The relative 
error of the period computed by this method is less at least by 4 order of 
magnitude in comparison with other numerical schemes. 
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Figure 1: TavgiNo,20) for the Surisl scheme (A^o < 3100), e 
Tth = 11, 65758528, T = 11, 88884005. 
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Figure 2: Tavg{No, 20) for the Surisl scheme (for very large Nq), e = 0.2, pQ = 1.95, 
Tth = 11,65758528, T = 11,88884005. 
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Figure 3: T/v for the Surisl scheme (for very large A^), e 
Tth = 11, 65758528, T = 11, 88884005. 
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Figure 5: Relative amplitude of the period oscillations (r) for po = 1, 8. Black cir- 
cles: symmetric projection, pluses: standard projection, white diamonds: discrete 
gradient, black squares: Surisl, stars: modified discrete gradient, black triangles: 
Suris2, black circles: leap-frog, dashed line: implicit midpoint. 
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Figure 6: Tn for the leap-frog scheme, e = 0.05, pp = 1.8, T = 9, 1254145545. 
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Figure 7: An for the leap-frog scheme, e = 0.05, pp = 1.8, T = 9, 1254145545. 
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Figure 8: Tn for the leap-frog scheme, e = 0.1, po = 0.05, T = 6,2815504224. 
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Figure 9: An for the leap-frog scheme, e = 0.1, po = 0.05, T = 6, 2815504224. 
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Figure 10: Tn for the Surisl scheme, e = 0.1, po = 0.05, T = 6, 297237955. 
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Figure 11: An for the Surisl scheme, e = 0.1, pp = 0.05, T = 6, 297237955. 
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Figure 12: Tavg = Ta„g(0, 100, 200) as a function of e for = 0, 1 {Tth = 
6,28711782). White squares: Surisl, black triangles: midpoint (Suris2 and dis- 
crete gradient methods yield practically the same results), black diamonds: mod- 
ified discrete gradient (very close to the theoretical exact values), white triangles: 
leap-frog, black circles: symmetric projection, white circles: standard projection. 




Fi gure 13: A^vg = ^aDg(0, 50) as a function of e for — 1,8 (^t/i — 2,239539). 
White triangles: leap-frog, black triangles: implicit midpoint, white squares: 
Surisl, black squares: Suris2, black diamonds: modified discrete gradient (dis- 
crete gradient method yields practically the same values), black circles; symmetric 
projection, white circles: standard projection (usually covered by black circles). 
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Figure 14: Modified discrete gradient method. Relative error as a function of e 
for po = 0, 02, Tth = 6, 283342395, T{e) = Tavg{0, 30). 
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Figure 15: Relative error as a function of e for po = 1.21 ("resonance value" 
for the leap-frog scheme), Tth = 7,01866131087, T{e) = fa„g(0, 100, 200). Black 
triangles: implicit midpoint, white triangles; leap-frog, black diamonds: modified 
discrete gradient. 
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Figure 16: for po = 2,000001, £ = 0,1. White triangles: leap-frog, black 
diamonds: modified discrete gradient, white diamonds: implicit midpoint. The 
period of the exact solution (continuous line): T^/j = 16,58809538, the average 
period given by the modified discrete gradient scheme: T = 16, 56380722. 




37 



Figure 17: Lpn for = 2, e = 0,2, round-off error A = 10~^^. White circles: 
standard projection, black circles: symmetric projection, white diamonds: discrete 
gradient, black diamonds: modified discrete gradient. 
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Figure 18: for pq = 2, e = 0, 00025, round-off error A = 10~^^. Black squares: 
Surisl, black circles: symmetric projection, white circles: standard projection, 
black diamonds: modified discrete gradient, white diamonds: discrete gradient. 
White squares (Suris2) and white triangles (leap-frog) are almost covered by black 
squares. 
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Table 1: Stability of the period. Minimal, maximal and average values of Tjv for po = 
1.95, g = 0.2 {Tth = 11.65758528). 



discretization 


maximal Tjy 


minimal Tjy 


average: T„„g (JV, 100, 200) 


N < 100 


JV ss 1.8 • 10® 


JV < 100 


JV fs 1.8 • 10® 


JV = 


JV = 1.8 ■ 10® 


leap-frog 


11.93166041 


11.93166040 


11.93164145 


11.93164140 


11.93165174 


11.93165162 


Surisl 


11.88885008 


11.88885008 


11.88883061 


11.88883061 


11.88884005 


11.88884001 


discrete gradient 


11.64698500 


11.64698540 


11.64697157 


11.64697190 


11.64697732 


11.64697764 



Table 2: Relative error of the amplitude {A = Aavg{0, 50)). 



PO 


leap-frog 


Surisl 


Suris2 


gradient 


mod. grad. 


projection 


sym. proj. 


midpoint 


e = 0.02 


0.05 


5.00E-05 


1.50B-04 


l.OOE-04 


-1.86E-08 


-1.87B-08 


-1.68E-08 


-1.71E-08 


-2.89E-08 


0.1 


6.00B-06 


1.50E-04 


9.99B-06 


-1.86E-08 


-1.84E-08 


-l.lOE-08 


-1.21E-08 


-6.01E-08 


0.3 


6.00B-06 


1.48E-04 


9.92B-06 


-1.74E-08 


-1.68E-08 


4.58E-08 


6.09B-08 


-3.95E-07 


0.6 


6.00B-06 


1.46E-04 


9.79B-06 


-1.66E-08 


-1.66B-08 


1.69E-08 


-8.69E-09 


-1.08E-06 


0.8 


6.02B-06 


1.39E-04 


9.47B-06 


-9.03B-09 


-8.37B-09 


-1.46E-07 


-2.06B-07 


-2.84E-06 


1.2 


6.13B-06 


1.26E-04 


8.86B-06 


-3.86B-09 


-3.88E-09 


-4.66E-09 


2.18B-09 


-7.00E-06 


1.6 


6.66B-06 


LOSE- 04 


8.24B-06 


2.71B-09 


2.68E-09 


5.42E-10 


1.33B-08 


-1.53E-05 


1.8 


6.73B-06 


1.02E-04 


8.48B-06 


4.07B-09 


3.96E-09 


4.56E-09 


4.10B-09 


-2.49E-05 


£ = 0.6 


0.05 


2.64B-02 


8.52E-02 


6.68B-02 


-6.34E-03 


-6.87B-03 


-3.16E-02 


-3.16E-02 


-6.35E-03 


0.1 


2.55E-()2 


8.52E-02 


5.57E-02 


-6.32E-03 


-6.8QE-03 


-3.05E-02 


-3.14E-02 


-6.34E-03 


0.3 


2.58E-02 


8.46E-02 


5.56E-02 


-6.07E-03 


-6.59E-03 


-2.55E-02 


-2.76E-02 


-6.27E-03 


0.6 


2.65B-02 


8.34E-02 


6.64B-02 


-6.72E-03 


-6.13E-03 


-1.44E-02 


-2.13E-02 


-6.30E-03 


0.8 


2.81E-02 


8.05E-02 


5.46E-02 


-4.58E-03 


-5.01E-03 


8.86E-03 


-5.54E-03 


-6.12E-03 


1.2 


3.07E-02 


7.42E-02 


5.32E-02 


-2.44E-03 


-2.69E-03 


3.40E-02 


1.95E-02 


-6.54E-03 


1.6 


3.84B-02 


6.52E-02 


6.19B-02 


1.80B-04 


1.46E-04 


9.31E-03 


1.76B-02 


-9.00E-03 


1.8 


4.76B-02 


6.14E-02 


6.46B-02 


1.22B-03 


1.31E-03 


5.56E-03 


6.70B-03 


-1.36E-02 
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Table 3: Relative error of the period (T = Tayg{0, 100, 200)). 



PO 


leap-frog 


Surisl 


SuriB2 


gradient 


mod. grad. 


projection 


sym. proj. 


midpoint 


e = 0.02 


0.02 


-1.67E-05 


8.33B-05 


3.33B-05 


3.33B-05 


-3.34B-09 


-1.66B-06 


-1.66B-06 


3.33E-05 


0.05 


-1.66E-05 


8.33B-05 


3.33B-05 


3.33B-05 


-2.08B-08 


-1.64B-05 


-1.66B-05 


3.33E-05 


0.1 


-1.66E-05 


8.32E-05 


3.33B-05 


3.32B-05 


-8.34B-08 


-1.56B-05 


-1.59B-05 


3.32E-05 


0.3 


-l.SOE-OF) 


8.18E-05 


,:i.3()E-05 


3.26E-05 


-7.52E-07 


-fi.74E-0fi 


-l.OlE-05 


3.24E-05 


0.5 


-1.45E-05 


7.92E-05 


3.23E-05 


3.12E-05 


-2.10B-06 


l.llE-05 


1.70E-06 


3.07E-05 


0.8 


-1.08E-05 


7.28E-05 


3.10B-05 


2.79B-05 


-5.45B-06 


5.58E-05 


3.16E-05 


2.63E-05 


1.0 


-6.99E-06 


6.71E-05 


3.01E-05 


2.47E-05 


-8.63E-06 


9.86E-05 


6.0,-)E-05 


2.20E-05 


1.2 


-1.48E-06 


6.05E-05 


2.95E-05 


2.07E-05 


-1.27E-05 


1.53E-04 


9.80E-05 


1.62E-05 


1.4 


6.86E-06 


5.37E-05 


3.03B-05 


1.56B-05 


-1.77B-05 


2.21E-04 


1.46E-04 


8.30E-06 


1.6 


2.12E-05 


4.92E-05 


3.52E-05 


9.31E-0fi 


-2.40E-05 


3.05E-04 


2.07E-04 


-3.63E-06 


1.8 


5.64E-05 


5.91E-05 


5.77E-05 


9.19E-07 


-3.24B-05 


4.08B-04 


2.87E-04 


-2.75E-05 


1.95 


2.17E-04 


1.90B-04 


2.03B-04 


-9.09B-06 


-4.24B-05 


4.99E-04 


3.72E-04 


-1.15E-04 


2.05 


-2.44E-04 


-2.78E-04 


-2.61B-04 


-1.14B-05 


-4.47B-06 


-3.47B-06 


1.57E-04 


1.14E-04 


2.2 


-9.25E-05 


-1.13E-04 


-1.03B-04 


-7.02B-06 


-4.04B-05 


-7.22B-05 


1.04E-04 


4.10E-05 


2.5 


-5.71E-05 


-6.97E-05 


-6.34B-05 


-4.20B-06 


-3.76B-05 


-1.08B-04 


6.77E-05 


2.54E-05 


3 


-4.45E-05 


-5.18E-05 


-4.81B-05 


-2.44B-06 


-3.58B-05 


-1.28B-04 


4.18E-05 


2.04E-05 


5 


-3.61E-05 


-3.83E-a5 


-3.72B-05 


-7.26B-07 


-3.41B-05 


-1.43B-04 


1.46E-05 


1.75E-05 


£ = 0.5 


0.02 


-1.06E-02 


5.07E-02 


2.05E-02 


2.05E-02 


-2.03B-06 


-1.06B-02 


-1.07B-02 


2.05E-02 


0.05 


-1.06E-02 


5.07B-02 


2.05B-02 


2.05B-02 


-1.25B-05 


-1.05B-02 


-1.06B-02 


2.05E-02 


0.1 


-1.06E-02 


5.06E-02 


2.0r)E-02 


2.04E-02 


-5.02E-05 


-9.86E-03 


-1.03E-02 


2.04E-02 


0.3 


-l.OlE-02 


4.97E-02 


2.03E-02 


2.01E-02 


-4.53B-04 


-3.29B-03 


-7.54B-03 


1.99E-02 


0.5 


-9.17E-03 


4.80B-02 


1.98B-02 


1.93B-02 


-1.27B-03 


l.OlE-02 


-1.69B-03 


1.89E-02 


0.8 


-6.71E-03 


4.40E-02 


1.89E-02 


1.73E-02 


-3.30E-03 


4.43E-02 


1.42E-02 


1.64E-02 


1 


-4.13E-03 


4.02E-02 


1.83E-02 


1.53E-02 


-5.25B-03 


7.85B-02 


3.12B-02 


1.38B-02 


1.2 


-4.05E-04 


3.58B-02 


1.79B-02 


1.29B-02 


-7.74B-03 


1.24E-01 


5.55E-02 


1.03E-02 


1.4 


5.31E-03 


3.11B-02 


1.83B-02 


9.82B-03 


-1.09B-02 


1.84E-01 


9.04E-02 


5.56E-03 


1.6 


2.40E-02 


3.74B-02 


3.08B-02 


8.57B-03 


-2.13B-02 


4.11E-01 


2.14E-01 


-1.91E-03 


1.8 


4.28E-02 


3.27B-02 


3.80B-02 


6.42B-04 


-2.03E-02 


3.1BE-01 


2.19E-01 


-1.56E-02 


1.95 


3.41E-01 


1.86B-01 


2.56B-01 


-5.72B-03 


-2.68B-02 


2.86E-01 


3.06E-01 


-5.94E-02 


2.05 


-1.17E-01 


-1.19E-01 


-1.19B-01 


-7.20B-03 


-2.83B-02 


-4.02B-02 


1.14E-01 


9.20E-02 


2.2 


-5.57E-02 


-5.95E-02 


-5.77E-02 


-4.45E-03 


-2.55E-02 


-4.42E-02 


8.08E-02 


2.70E-02 


2.5 


-3.68E-02 


-3.87E-02 


-3.77B-02 


-2.68B-03 


-2.37B-02 


-5.03B-02 


5.49E-02 


1.64E-02 


3 


-2.96E-02 


-2.96E-02 


-2.95B-02 


-1.57B-03 


-2.26E-02 


-6.61B-02 


3.34E-02 


1.34E-02 


5 


-2.68E-02 


-2.31E-02 


-2.50B-02 


-5.04B-04 


-2.14B-02 


-5.54B-02 


4.32E-03 


1.29E-02 
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Table 4: Relative error of the period in the neighbourhood of the separatrix. The blank 
space with a dot means that the qualitative behaviour of the discretization is wrong. 



PO - 2.0 


leap-frog 


Surisl 


Suris2 


gradient 


mod. grad. 


projection 


sym. proj. 


midpoint 


E = 0.02 


-l.OE-02 


8.96E-04 


8.50B-04 


8.73E-04 


-1.50E-05 


-4.83E-()5 


5.19E-()4 


4.08E-04 


-4.57E-04 


-l.OE-03 


7.12E-03 


7.06E-03 


7.09B-03 


-1.95E-05 


-6.29B-05 


5.14E-()4 


4.26E-04 


-3.40E-03 


-l.OE-04 


9.17E-02 


9.16B-02 


9.16E-02 


-2.22E-05 


-5.56E-05 


5.05E-04 


4.34E-04 


-2.40E-02 


-l.OE-05 








-2.43E-05 


-5.58E-05 


4.99E-04 


4.40E-04 


-1.03E-01 


-l.OE-06 








-2.80E-05 


-6.69B-05 


4.94E-()4 


4.43E-()4 


-2.13E-()1 


-l.OE-07 








-7.33E-05 


-2.09B-05 


4.91E-()4 


4.46E-()4 


-3.08E-01 


-l.OE-08 








1.38B-04 


1.15E-()4 


4.88E-()4 


4.48E-()4 


-3.83E-()1 


-l.OE-09 








-1.61E-03 


1.18E-()3 


4.86E-()4 


4.5()E-()4 


-4.43E-()1 


l.OB-08 


-4.16E-01 


-4.15E-01 


-4.15E-01 


-5.16E-05 


-4.23B-()6 


2.82E-04 


3.32E-()4 




l.OE-07 


-3.44B-01 


-3.44E-01 


-3.44E-01 


-1.59E-05 


-6.26B-05 


2.60E-()4 


3.15E-04 




l.()E-()6 


-2.54E-()1 


-2.54E-01 


-2.54E-01 


-2.90E-05 


-6.44E-05 


2.31E-04 


2.94E-04 




l.()E-()4 


-4.2fiE-()2 


-4.27E-()2 


-4.27E-()2 


-2,22E-()5 




8.9:iE-(),'') 


1.7fiE-04 


3.38E-02 


l.UE-u:! 


-(j.(jhE-u:i 


-(>.7:!i;-o:! 


-(i.7()E-l):! 


-i,!)(.E-0.-. 


-,").2!jE-iJ,") 


E()2E-U4 


2.Ij!)E-U4 


:!.4!)E-03 


l.OE-01 


-1.46E-04 


-1.74E-04 


-1.60E-04 


-9.25E-0(i 


-4.26B-06 


-3.91B-06 


1.31B-04 


6.62B-06 


£ = 0.5 


-l.OE-02 








-9.51E-03 


-3.06B-02 


2.23E-()1 


3.31E-()1 


-1.54E-()1 


-l.OE-03 








-1.24E-02 


-3.36E-02 


1.59E-()1 


3.31E-()1 


-3.21E-()1 


-l.OE-04 








-1.41E-02 


-3.53E-02 


1.19E-()1 


3.2fiE-()l 


-4.48E-()1 


-l.OE-05 








-1.52E-02 


-3.65B-()2 


9.()9E-()2 


3.22E-()1 


-5.36E-()1 


-l.OE-06 








-1.61E-02 


-3.73B-()2 


7.16E-()2 


3.19E-()1 


-6.()1E-()1 


-l.OE-07 








-1.67E-02 


-3.80B-02 


5.87E-02 


3.17E-01 


-6.49E-01 


-l.OE-08 








-1.73E-02 


-3.86B-02 


4.46B-02 


3.16B-01 


-6.88B-01 


-l.OE-09 








-1.73E-02 


-3.86B-02 


3.55E-02 


3.14E-01 


-7.18E-01 


l.OE-08 


-7.24B-01 


-7.22E-01 


-7.23E-01 


-1.72E-02 


-3.85B-()2 


-5.68E-()2 


2.42E-()1 
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